Nucleation and Growth of the Superconducting Phase in the Presence of a Current 
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We study the localized stationary solutions of the one-dimensional time-dependent Ginzburg- 
Landau equations in the presence of a current. These threshold perturbations separate undercritical 
perturbations which return to the normal phase from overcritical perturbations which lead to the 
superconducting phase. Careful numerical work in the small-current limit shows that the amplitude 
of these solutions is exponentially small in the current; we provide an approximate analysis which 
captures this behavior. As the current is increased toward the stall current J*, the width of these 
solutions diverges resulting in widely separated normal-superconducting interfaces. We map out nu- 
merically the dependence of J* on u (a parameter characterizing the material) and use asymptotic 
analysis to derive the behaviors for large u (J* ~ m"^''^) and small u (J ^ Jc, the critical deparing 
current), which agree with the numerical work in these regimes. For currents other than J* the 
interface moves, and in this case we study the interface velocity as a function of u and J. We find 
that the velocities are bounded both as J ~-+ and as J —> Jc, contrary to previous claims. 

PACS numbers: 74.40. -Hk, 74.60. Jg, 74.76.-w, 64.60. Qb 



I. INTRODUCTION 

When a superconductor is placed in a magnetic field 
equal to its critical field He, the normal and supercon- 
ducting phases can coexist in a state of equilibrium with 
the two phases separated by normal-superconducting 
(NS) interfaces. The dynamics of such interfaces is im- 
portant for various nonequilibrium phenomena. For in- 
stance, if the applied magnetic field is quenched be- 
low He, these interfaces move through the sample, ex- 
pellin|S the magnetic flux so as to establish the Meissner 
phasad^Ei. Just as superconductivity can be destroyed 
by applying a magnetic field exceeding He, it can also 
be destroyed by applying a current exceeding the crit- 
ical depairing current Jc- Thus by analogy with the 
magnetic field case, one might expect a special value of 
the applied current J* < Jc at which the superconduct- 
ing and normal [phases coexist, separated by a station- 
ary NS interfaceliy. In contrast to the magnetic field in- 
duced NS interfaces, these current-induced NS interfaces 
are intrinsically nonequilibrium entities, and their struc- 
ture depends upon the dynamics of the order parameter 
and magnetic field. The evolution and dynamics of these 
nonequilibrium interfaces is the subject of this paper. 

The current-induced NS interfaces arise in several con- 
texts. First, they are known to be important in under- 
standing the dynamics of the "resistive statej!L.in super- 
conducting wires and films (for a review seeEJ), and in 
determining the global stability of the normaLaiid super- 
conducting phasjcs in the presence of a currenttHI. Second, 
Aranson et aZJl3 have recently studied the nucleation of 
the normal phase in thin type-II superconducting strips 
in the presence of both a magnetic field and a transport 
current. They found that a sufficient current produced 
large normal droplets containing multiple fiux quanta. 



Without a current one finds stationary, singly quantized 
vortices, with a larger amount of NS interface per flux 
quantum than a multiply quantized droplet. They con- 
clude that the current produces an effective surface ten- 
sion for the NS interface which is positive, stabilizing 
the interface and producing larger droplets with smaller 
surface area. Motivated in part by its role in this phe- 
nomenon we wanted to re-examine the nonequilibrium 
stabilizing effects of current. 

Even when the superconducting phase is ostensibly the 
equilibrium phase, a current makes the normal phase 
metastable, i.e., linearly stable to infinitesimal supercon- 
ducting perturbations. A localized superconducting per- 
turbation of finite amplitude, on the other hand, has one 
of two fates: (1) its amplitude may ultimately shrink 
to zero restoring the normal phase (undercritical) or (2) 
it may grow eventually establishing the superconduct- 
ing state (overcritical). Separating these two possibil- 
ities are the critical nuclei or threshold perturbations, 
for present purposes stationary solutions of the time- 
dependent Ginzburg-Landau (TDGL) equations local- 
ized around the normal state. As one raises the current, 
the amplitude of the threshold solution grows, implying 
that the normal phase becomes increasingly stable. 

At very low currents, the widths of the critical nuclei 
shrink as the current is increased, but eventually this 
trend is reversed and the width grows as the current is 
increased further. In fact, as the current-approaches a 
particular value, J* (the "stall current"ll3), the width 
diverges resulting in two well-separated, stationary NS 
interfaces. The interfape solutions have been studied nu- 
merically by Likharevlij, who found that the interfaces 
were stationary at J* « 0.335 for u — 5.79, where u 
characterizes the material and is 5.79 for nonmagnetic 
impuritiesEj. They were also studied by Kramer and 
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Baratoffl, who found J* « 0.291 for. u = 12 (corre- 
sponding to paramagnetic impuritiescJ) . However, we 
know of no systematic study of the dependence of J* 
upon u. In this work we remedy this situation by us- 
ing a combination of numerical methods j-and analysis 
including matched asymptotic expansionsE3. We show 
that J* for large u in contrast to a previous 

conjectureEHI, and we find how J* approaches Jc in the 
small-u limit. 

At currents other than J*, the interfaces move with 
a constant velocity with the superconducting phase in- 
vading the normal phase for J < J* and vice versa for 
J* < J < Jc- At currents close to J*, the-interface ve- 
locity is proportional to {J ~ J*). Likharevll3 defined the 
constant of proportionality, rj — {dc/dJ)~^\j=j* , where 
c is the interface speed; he found 77 w 0.7 for u = 5.79. 
In the extreme limits, J — s- and J ^ Jc Likharev pre- 
dicted that the speed c diverges. We find c to be bounded 
in both cases and provide an analytic expression for it as 
J-^ 0. 

The results of this work are summarized in Table |. 
The rest of the paper is organized as follows. After briefly 
reviewing the TDGL equations and the approximations 
used in this work (section |l), we study the critical nuclei 
focu sing on their size and shape in the limit J — > (sec- 



tion III). We then move on to consider the stationary 
interface solutions; in particular we map out the depen- 
dence of the stall current J* on u and supplement the 
numerical work with analysis of the u — > cx) and u ^ 
limits (section IV). Next, we examine moving interfaces 
first in the linear response regime and then in the limits 
J ^ and J ^ Jc (section ^ . The appendix contains a 
calculation of the amplitudes of the critical nuclei in the 
J ^ limit. 

TABLE I. Summary of the primary results. 



I. Critical nucleus 

Small- J width 
Small-J amplitude 

II. Stall current J* 
Large-w 

Small-M 

III. Kinetic coeff. "q 
Large-u 
Small-u 

IV. Interface speed 
J ^ 



W 



' exp{— A/uJ} 



J* = 0.584491 u-^^* 
(1- V8)^/^ 



J 



Jc 



(1 - u/24) 



3/2 



= 0.797 u'^^'^ 



* 2/u 

„ -0.85 



-0) 

{u —> 00) 
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Sec. 

Sec. 
Sec. 

Sec. 
Sec. 
Sec. 
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II. THE TDGL EQUATIONS 

The starting point for our study is the set of TDGL 
equations for the order parameter ip, the scalar potential 
$, and the vector potential A: 



nj[dt + 



2m 

47r , 



ie*A 

he 



V X V X A = — (J„ + J,), 

c 



(2.1) 



(2.2) 



where the normal current J„ and the supercurrent Js are 
given by 



J„ = (t(") (-atA/c-V$), 



(2.3a) 



Tip* p*2 

J. = ^(^*V^-VV^*)-— IV'PA, (2.3b) 

and where 7 (which is assumed to be real) is a di- 
mensionless quantity characterizing the relaxation time 
of the order parameter, cr(") is the normal conductiv- 
ity, and a — ao(l — T/Tcq). From these parameters 
we can form two important length scales, the coher- 
ence length ^ ~ fi/{2m\a\y^^ and the penetration depth 
A = [m6cV47r(e*)2|a|]i/2. 

These equations assume relaxational dynamics for the 
order parameter as well as a two-fluid description for the 
current. With somewhat restrictive assumptions, thpy 
can be derived from the microscopic BCS theoryolij. 
Further simplification is possible in the limit of a thin, 
narrow film, that is, when the thickness is less than the 
coherence length, d < ^, aud the width is less than the 
effective penetration depthllZl, w ^ X^/d. In this case 
the current carried by the film or wire is small, and we 
needn't worry about the fields it produces. Therefore, 
Eq. (2.2) may be dropped, and we need only specify the 
total current J = J„ -t- (subject to V • J = ), along 
with the order-parameter dynamics, Eq. ( |2.l| ). This 
approximation is commonly used for superconducting 
wireaij and cauj-be justified mathematically for super- 
conducting filmalS. In addition, we will be considering 
processes in the absence of an applied magnetic field, so 
that we may set A = 0. With these simplifications, we 
can now rewrite the equations in terms of dimensionless 
(primed) quantities. 



$ = 



;ie*|o| 



J _ / 2 e-|a|^^" r/ 
■J — \l ^ ft J , 



(2.4) 



which leads to 
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u{dt' + ifi'yy = (V'2 + 1 - 
J' ^ Im(7/;'*W) - Vy, 



V- J' 



0. 



(2.5a) 
(2.5b) 



Note Aat length is measured in units of coherence 
lengtho. We wiU drop the primes hereafter. The only pa- 
rameters remaining in the problem are the scaled current 
J and a dimensionless material parameter u = t^/tj, 
where = h,'y/\a\ is the order-parameter relaxation time 
and Tj — (T^")TO5/e*^|a| is the current relaxation time. 
We will treat m as a phenomenological parameter and 
study the nucleation and growth process as a function of 
u. The microscopic derivations of the TDGL equations 
predict that u — 5.79 (nonmagnetic impurities)lij, and 
u = 12 (paramagnetic impurities )ll3, but sraall u is also 
useful for modeling gapped superconductorsEJ. 



III. NUCLEATION OF THE 
SUPERCONDUCTING PHASE FROM THE 
NORMAL PHASE 

In the presence of an applied current the normal phase 
in a wire is linearly stable with respect to supccmiiduct- 
ing perturbations for any value of the currentEiJO. The 
reason for this stability is that any quiescent supercon- 
ducting fluctuation will be accelerated by the electric 
field, its velocity eventually exceeding the critical depair- 
ing velocity, resulting in the decay of the fluctuation. 
The growth of the superconducting phase therefore re- 
quires a nucleus of sufhcicnt size that will locally screen 
the electric field and allow the superconducting phase to 
continue growing; smaller nuclei will simply decay back 
to the normal phase. The amplitude of the "critical" 
nucleus should decrease as the current approaches zero, 
reaching zero only at J = 0. We expect the critical nuclei 
to be unstable, stationary (but nonequilibrium) solutions 
of the TDGL equations, which asymptotically approach 
the normal solution as x — > ±oo. These "bump" so- 
lutions of the TDGL equations are the subject of this 
section. We include here an extensive numerical study of 
the amplitudes and widths of the critical nuclei, as well 
as some analytical estimates for these quantities. 



A. Numerical results 

Let us start by discussing the numerical work on the 
critical nuclei. For the analytic work, we often find it 
convenient to use the amplitude and phase variables, i.e. 
"0 = /e*^; but they are ill-suited for the numerical work, 
since the calculation of the phase becomes difficult when 
the amplitude is small. Following LikharevlllJ we use in- 
stead Tp = R+il, with R and / real, and in one dimension 
Eqs. ~ 



(2.5) become 

uRt = Rxx + u^il 
ult = Ixx - UfJ.R - 

J = RIt — IRr 



R-{R^ 
I-{R^A 

- IJ^x- 



(3.1a) 
(3.1b) 
(3.1c) 



Since the nuclei are unstable stationary states, they 
are investigated only by time-independent means. Such 
solutions require a particular gauge choice — in this case 
/i(a;) = where il^i^) has its maximum amplitude; the; 
are then sought using a relaxation algorithnJHj. Figure 
shows a typical bump's amplitude, / = \/R^ + P, the 
associated superfluid velocity q = {RIx — IRx)/ f'^ and 
the electric field E{x) = —Hx{x)- The figure shows only 
half of the solution; f{x), q{x) and E{x) are even about 
a; = 0. In Fig. ||we plot the bump's maximum amplitude, 
-00, as a function of J; it grows as the current rises, in- 
dicating the increasing stability of the npcmal phase. In 
the data presented by Watts- Tobin et al. cfl appears to 
vary linearly with J for small J. However, in our numer- 
ical calculations at very small currents the dependence 
deviates from linearity (see the inset of Fig. H), and ipQ 
drops rapidly to zero as J — > 0, cons iatej it with the expo- 
nential behavior suggested in Refs.Elrnd. More precisely 
our small-J data (0.008 < J < 0.015) at u = 5.79 is fit 
by 



il^oiJ) = Bexp(-A/wJ), 



(3.2) 



with A = 0.042 and B = 0.19. A somewhat similar 
dependence (with A = 2/3) was suggested by Ivlev et 
a/.tifE3; they were considering a distinct quantity but 
one also related to critical fluctuations about the normal 
phase (see the Appendix for more details). 
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FIG. 1. The bump's amplitude, /(x), its superfluid veloc- 
ity, q{x), and the electric field, E{x), for u = 5.79 and J — 0.2. 
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FIG. 2. The maximum amplitude of the bumps tpo as a 
function of J for u = 5.79. The inset shows the exponential 
dependence of the small-J data, see Eq. (3.2). 



The width of the bump diverges in the small-J limit 
like (mJ)~^/^, as can be seen from the analysis below. So 
as J is increased from zero, the width initially shrinks, 
but eventually the width begins to grow again, diverg- 
ing as the current approaches the stall current J*. In 
this limit the bump transforms into two well separated 
interfaces (see Fig. 0). 




5 10 15 

Distance from center of bump, x 

FIG. 3. The bump profiles for J* - J = 10"^ (sohd), 
* - J = lO''^ (dashed) and J* - J = 10"^ (dotted) at 
= 5.79. 



where we have dropped the term ipt and selected the 
gauge /Lt(0) = 0. We saw in Fig. ^ that ipQ becomes very 
tiny in the small-J limit, thus the nonlinear terms can 
be neglected, leading to 



Tpxx + iuJxtp + ^p = 0, 



(3.4) 



a complex version of the Airy equation. Applying the 
WKB method results in the approximate solution 



[1 + 
xexp 

xexp 



(uJx)- 
2 

3 uJ 



-1/8 

fi , / T n2i3/4 3a 

[1 -I- {uJxYl COS — 



3 uJ 



[l + (uJx) 



21 3/4 



3a 



(3.5) 



where a ~ tan^^(wJa;). The numerical data agrees quite 
well with this predicted shape in the small-J limit. For 
small X the expression can be approximated by 



■0 exp [i(l - uJ l^)x ~ uJx"^ /A\ 



(3.6) 



We see here that the width of the bump varies like 
(mJ)""'^/^ in this limit and that the superfluid velocity 
q K, {1 ~ uJ/A). For large x, on the other hand, where 
a ~ 7r/2, the expression becomes 



exp 



i^r(i-z) 



(3.7) 



as one expects for the Airy function. Note that deep in 
the tail of the solution, we see a different length scale 
^Airy ^ {uJ)~'^/^ arising. 

Since the above analysis is of a linear equation, it can 
not determine the amplitude of the nucleus; for this pur- 
pose the nonlinearities must be considered. In the ap- 
pendix we outline an ad hoc calculation of the small-J 
limit of the bump amplitude. We take a ■(/' of an un- 
known amplitude but of a fixed shape inspired by the 
above analysis and assume that it is a stationary solu- 
tion of the full TDGL. We then determine its amplitude 
self-consistently. The resulting amplitude is 



1^1 



2J 



1/4 



-1/2 



exp 



16 



81 uJ 



(3.8) 



The factor A = 16/81 is within a few percent of that 
extracted from the numerical data. 



B. Analysis in the J ^ limit 



IV. STATIONARY INTERFACES 



The equations for nuclei centered at the origin are 

'4'xx — iufiip + il^ — = 0; (3.3a) 

li = —Jx+ / ln\{'ip* ij^x') dx , (3.3b) 
Jo 



As the current is raised, the width of the critical nu- 
cleus grows and ultimately diverges as the stall current 
is reached, resulting in well separated, stationary inter- 
faces. These interface solutions will be the subject of the 
rest of this work. 
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A. Numerical methods and results 



Let us first discuss the numerical work on the interface 
solutions. For given values of u and J we evolved the 
TDGL equations from an initial guess which is purely 
superconducting on the left, iIj{x) = /oo e**°°^ and 
fJ'x{x) = 0, and purely normal on the right, ipix) = 
and ^x{x) — —J. The values /oo and q^o are related to 
the applied current through 



2 

00 ' 



(4.1a) 
(4.1b) 



Stability requires^ taking the larger positive root of the 
former equations which places the following bounds on 
J, /oo and goo: 



< J < Jc = ^4/27 w 0.3849, 

1 > /oo > v/273 « 0.8165, 
< goo < 7173 « 0.5774. 



(4.2a) 
(4.2b) 
(4.2c) 



We employed several schemes to integrate the equa- 
tions in time includipg both explicit (Euler) and implicit 
(Crank-Nicholson)Ej. 

Initially the front moves and changes shape but even- 
tually it reaches a steady state in which the interface 
moves at a constant velocity without further deforma- 
tion. By the time-dependent means we found locally sta- 
ble, constant-velocity solutions for currents less than J^. 
To examine these solutions more accurately, we adopted 
a time-independent method. First, we transformed coor- 
dinates to a moving frame, x' = x — ct; next, we chose 
fJ- = cgoo as x — !■ — 00 which allows for a truly time- 
independent solution (i.e. one with both amplitude and 
phase time-independent). Then wc searched for station- 
ary solutions using a relaxation algorithmE3 where {u, J) 
are input parameters and c is treated as an eigenvalue. 
This approach requires an additional boundary condition 
to fix translational invariance; we elected to fix fj, on the 
rightmost site. To find the stall currents J* we can set 
c = and take J or u as the eigenvalue. 

Figures ^ and |^ show the order-parameter amplitude / 
and the electric-field distribution E — —ji^ of the stalled 
interface determined for u = 500 and u ~ 1.04 respec- 
tively. Note that while / is very flat in the supercon- 
ducting region, the real and imaginary parts, R and /, 
oscillate with a wavelength 27r/goo- Because of this addi- 
tional length scale inherent in R and /, there is little to 
be gained from varying the mesh size. In fact, this length 
scale is compressed as we move to the right, and we are 
only saved from the difficulties of handling rapidly oscil- 
lating functions by the fact that the amplitudes decay so 
quickly. 
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FIG. 4. The stationary NS interface solution when u — 500 
for which the stall current J* = 0.12252. Shown here are 
the numerically determined f{x) and E (x), as well as the 
Langer-Ambegaokar (LA) solution (Eq. (4.10), the solution 



with no electric field) corresponding to the same current. 
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FIG. 5. The stationary NS interface solution when 
u = 1.04 for which the stall current J* = 0.3836. Shown here 
are the numerically determined f{x) and E{x), as well as the 
function fo{x), the u ^ profile, derived from Eq. (4.27) 



In the large-u case (see Fig. ^), E{x) remains flat 
throughout most of the space; it changes abruptly from 
one constant to another only after f{x) has become small. 
The variations of f{x) are more gradual; however, the 
greatest changes in occur in that same small area. 
This region of rapid change is known as a boundary layer, 
it marks where the current suddenly changes from super- 
conducting to normal, i.e., the position of the NS inter- 
face. As u increases, the longer length scale over which 
/ varies on the superconducting side remains essentially 
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fixed, while the boundary-layer thickness shrinks to zero. 
In the opposite limit, the small-u case (see Fig. ||), f{x) 
and E{x) appear to vary together even in the supercon- 
ducting region; moreover, the length scale over which 
they vary grows as u is decreased. We will postpone 
providing more of the numerical results on the interfaces 
until some of the analytic arguments are available for 
comparison. 



B. Asymptotic analysis of the interface solutions: 
preliminaries 

Before addressing the large-w and small-u limits sepa- 
rately, let us put the TDGL equations into a form con- 
venient for analysis and derive expressions for the length 
scales deep in the superconducting and normal regions. 
The disparity of these length scales in the large-u limit 
will motivate the boundary-layer analysis in that regime; 
while an inequality they satisfy will lead to the conclusion 
that J* — *■ Jc in the small-u limit. 

We make the substitution tp = /e*^, which yields 



uft = - fie.)^ + f- 
u{0t + ^l)f = 2fj, + je. 



(4.3a) 
(4.3b) 
(4.3c) 



Next we restrict our attention to stationary solutions. 
Note that only spatial derivatives of 9 appear now, al- 
lowing us to work with the superfluid velocity q — 9x 
instead of 9. The equations become 



'/ + / 



J XX 



fqx, 



0, 



^J'x 



(4.4a) 
(4.4b) 
(4.4c) 



where J* replaces J as these e quati ons apply to the stall 
situation. Next multiply Eq. ( 4.4b| ) by / and note that 
the right hand side is now {f'^q)x which we can express in 
terms of ^ by differentiating Eq. ( 4.4c ); these steps lead 
to 

IJ-xx = upn. (4.5) 
Now let us assume the following asymptotic forms as 



lim f{x)^ f^- ue'^i + 



lim q{x) = Qa, 

c — * — oo 

lim fJ.{x) = /io 



qi e 



x/\q 



x/\fi 



(4.6a) 
(4.6b) 
(4.6c) 



Subs titut ing these expressions into Eqs. (4.4a), (4.4c) 
and (|4.5| ) and recalling the boundary conditions yields 



(2fl - Xj^) he^/^f - - 0, (4.7a) 

-2/o,goo/ie"/^^ + flqie^^^- + A" Vie"/^- = 0, (4.7b) 



- nfl = 0. 



(4.7c) 



Eq. ( 4.7c ) provides an expression for A^, the electric-field 
screening length. Since foo is always of 0(1), we see that 
A^ shrinks as u — > cx) and diverges as u ^ 0, which is 
consistent with the behavior seen in Figs.|4|and|,_^ 

Mo re than one decay length appears in Eqs. ( 4.7a ) and 
( 4.7b| ). If they are not equal, the term with the shorter 
length is exponentially small compared to the other(s) 
and will not contribut e to t he a: — > — oo limit. Since none 
of the terms in Eq. ( 4.7b| ) can equal zero individually, 
we conclude that the longer two of A/, \q and A^ must 
be equal. Next, because the term multiplying e^/'^i in 
Eq. (4.7a) cannot equal zero on its own, we determine 
that \q < A/, making Ay one of the longer lengths. Fi- 
nally, if we assume that A/ A,, 



> Xq we find that 



A/ = and A^ = u^^^^f^^ and reach a contra- 

diction (except at u — 2). Thus provided the original 
assumption of an exponential approach is valid, we con- 
clude that 



A 



/ 



Aq > A^. 



(4.8) 



This equality of A/ and Xq is reasonable given that both 
/ and q are related to the complex order parameter ip. 
Also, having A/ > A^ is consistent with the large-u data 
seen in Fig. |[ If A^, 7^ A/ then 



A 



/ 



6/o 



A 



LA- 



(4.9) 



We identify this length scale as Xla since it coincides 
with that occurring in the solution of Eqs. (4.4) without 
any electric field (/i(a;) — 0), 




x , (4.10) 



fix)^fl-i'3fl-2) sech^ 



which was found by Langer and Ambegaokar@ in their 

study of phase slippa ge. The asymptotic form of 

Eq. ( [4.10| ) looks hke Eq. ( |I.6aD with A/ given by Eq. (p|). 



As a matter of fact because Ay ^ A^ in the large-u limit, 
the profile of f{x) is only imperceptibly different from 
the Langer- Ambegaokar (LA) solution in the supercon- 
ducting region and deviates from it only in the boundary 
layer, as is shown in Fig. IJ. 

Recall that A^ diverges as u — > 0; the inequality 
A/ > Ap implies that A/ must diverge as fast or faster 
in this limit. This scenario is consistent with the small-u 
data shown in Fig. |^ in whi ch f{x) and E(x) vary on 
long length scales. Eq. ( [4.9|) suggests that a diverging 



Xf implies that 



■\/2/3 and in turn that J Jc 
as u — !■ 0, which is also consistent with what is found 
numerically. 

In the other asymptotic limit, deep in the normal 
regime, tp is very small and hence the nonlinear terms 
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in Eqs. (21) can be dropped as was done for the bumps 
in the smaU-J hmit. The result is a complex Airy 
equatio n, th e asymptotic analysis of which was supplied 
in Eq. (3/?), where we saw the length scale XAiry ^ 



-> oo 
The 



(mJ*)~^/^. Somewhat like \Airy shrinks as u 
and expands as m — > but with different powers of u. 
presence of the disparate length scales, A/, and \A^ry^ 
in the large-w limit motivates the use of the boundary- 
layer analysis that comes next. We will see that XAiry 
scales in the same way as the boundary-layer thickness. 



we need to ensure that the solutions in the boundary layer 
match onto the solutions in the superconducting and nor- 
mal regions. By expanding the superconducting solution 
near the interface, we see that f{x) ~ —x/\/2 as the 
boundary layer is approached; matching to the boundary 
layer requires /j; '--^ 1, so that a — ^. In the normal re- 
gion, ^ « —J*x, so that matching to the boundary layer 
requires fix J* , and P = j + S. Solving this set of equa- 
tions, we conclude that a — j ~ 5 ~ 1/A and (3 = 1/2, 
i.e., the stall current J* ^ u^^l^ for large u. Note that 
J* as u — !■ cx), so tha t we were justified in drop- 
ping J*V/3 from Eq. Substituting J* - u-^l^ 



C. Asymptotic behavior of the stall current as it ^ cxj into XAiry gives A^ 



iiry 



We have already seen in Fig. ^ that the large-w pro- 
file can be divided into two regions — one slowly varying, 
one rapidly varying, also known as the outer and inner 
regions, respectively. Furthermore, it has been suggested 
that the ratio of the length scales characterizing these re- 
gions decreases as u — + oo. These features make the prob- 
lem ideally suited for boundary-layer analysis, in which 
one identifies the terms that dominate the differential 
equation in each region, analyzes the reduced equations 
consisting of dominant terms and then matches the be- 
havior in some intermediate region. 

We sta rt by eliminating the superfluid velocity q from 
Eqs. (4^), resulting in 



fxx 
f^xx 



{J* 



+ fixfr^ + f-f = o, 



uffl = 0. 



(4.11a) 
(4.11b) 



Let us consider first the slowly varying, superconducting 
region. We saw in the preliminary analysis that for large 
M, fJ-{x) is exponentially small, so we drop it. Next, let us 
assume that J* is small and drop it; we can verify in the 
end that this is self-consistent. The reduced equation is 



fxx ^" / / 



0, 



(4.12) 



with solution f{x) = — tanh(.T/\/2). 

Moving in from the left toward the interface (into the 
boundary- layer region), / becomes small, and the sec- 
ond term in Eq. (4.11a) which was subdominant becomes 
dominant. In this inner region / is small but rapidly 
varying, thus the dominant terms are 



[J* +flxf 

p 



(4.13) 



along with Eq. ( 4.11b ). Having identified the dominant 
terms, now we must make certain they balance. We as- 
sume that in the boundary layer, all the quantities scale 
as powers of u: 



,-/3 



J* 



-5 



(4.14) 



Balancing terms in Eq. ( 4.13 ), we find 2q; = 7 + J, while 
balancing terms in Eq. (4.11b) yields 2{a + 5) = 1. Next, 



indicating that it may be 
identified as the boundary-layer thickness. 

In order to determine the coefhcient of the u~^^^ in 
the stall current we need to solve the boundary layer 
problem. Let us rescale in the way suggested above: 



/ = u-i/^F, fi = u-^'^M{X), 



(4.15) 



Substituting these rescaled variables into Eqs. (4.11a 
and ( l.lllj) , and then expanding F, M and J in pow- 



ers of u we obtain at the lowest order 



o,xx 



{Jo + Mo,x) 
Mq,xx - F^Mo = 0, 



= 0, 



(4.16a) 
(4.16b) 



with the boundary conditions (from the outer regions) 

Fo,x(-oo) = -l/V2, Mo(-oo)=0, (4.17) 
Fq{+(x) =0, Mq,x{+(x) =-Jq. (4.18) 

(As before we need an extra boundary condition to fix 
the translational invariance.) For a n arbitrary Jo the 



solutions of Eqs. ( 4.16a ) and ( 4.16b ) cannot satisfy the 
boundary conditions; Jg must be tuned to a particular 
value before all of the boundary conditions are satisfied, 
leading to a nonlinear eigenvalue problem for Jq. We 
have solved this eigenvalue problem numerically and find 
that Jo = 0.584491. Therefore, to leading order we have 
for the stall current 



J* w 0.584491 u"^/'*. 



(4.19) 



This prediction agrees well with the numerical results 
and disagrees with Likharev's conjecture of a u~^^^ 
dependencelij, as can be seen in Fig. I and in Table || It 
is in principle possible to carry out this procedure to suc- 
cessively higher orders, but the equations become cum- 
bersome. Instead we have simply opted to fit our numer- 
ical data to a form inspired by the asymptotic analysis, 

J* = 0.584491 - 0.117461 u"^/^ - 0.12498 

-1-0.163043 u"^/'' + O (u-^/^) . (4.20) 
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D. Asymptotic behavior of the stall current as « ^ 



0.8 - 
0.6 - 

0.4 - 



0.04 - 



10° 



Numerical Results 
0.584491 u"" 



10' 



10=^ 



10^ 



10" 



10= 



FIG. 6. A log-log plot of the stall current J* vs. u. The 
solid line shows the numerically determined J*'s as a func- 
tion of u and the dotted line is 0.584491 u~^^^ (the large-it 
behavior predicted by matched asymptotic analysis). 



Now let us examine the opposite limit of u — > 0. In this 
case the eleGtric-field screening length becomes long, and 
Ivlev et aZ.E3 have proposed that this makes the small- 
it limit useful for modeling gapped superconductors. As 
already suggested the inequality of length scales, Xf > Xf_i 
implies that J* Jc- We will begin our small-M analysis 
by putting this result on firmer ground and extracting as 
a byproduct the u — )■ limit of the interface profile. 

The rescaled equations. Recall that deep in the super- 



conducting region 
that we rescale distance: 



,-1/2 



This observation suggests 
furthermore, to 



ensure that the normal current {—^x) scales in the same 
way as the total current we rescale /i = u^^^^ji as well. 
These rescalings yield 



J* = - Ai, 



0, 



(4.21a) 
(4.21b) 
(4.21c) 



TABLE 11. Representative numerical results for the stall 
current J* and kinetic coefficient rj. 



u 


J* 








1 


0.3838 


0.3838 


0.01871 


0.01871 


5 


0.3407 


0.5094 


0.6400 


0.1914 


10 


0.3013 


0.5359 


1.573 


0.2797 


50 


0.2127 


0.5655 


8.258 


0.4315 


100 


0.1807 


0.5715 


15.59 


0.4931 


500 


0.1224 


0.5788 


62.51 


0.5875 


1000 


0.1033 


0.5807 


111.3 


0.6259 


5000 


0.0693 


0.5828 


407.9 


0.6847 


10000 


0.0583 


0.5833 


708.4 


0.7084 


50000 


0.0391 


0.5840 


2487 


0.7440 



placing the small parameter u in front of f^x- If we ex- 
pand these functions as series in powers of u 



f = .fo + uh + ..., 

q = qo +uqi + . . . , 
/i = /to + Ufli + 
J* = Jq + tt Ji + . . . , 

then we find at the lowest order 
-fofo + fo- /o' = 0, 

Ao/o = 2^0/0,2; + Io<10,xi 



(4.22a) 
(4.22b) 
(4.22c) 
(4.22d) 



(4.23a) 
(4.23b) 
(4.23c) 



The solution of Eq. ( 4.23a ) is either fo — (the normal 
phase) or /o = (1 — fZo)^^^ i^^^ superconducting phase). 
Let us focus on the superconducting solutions. By elim- 
inating qo, we obtain the first order equations 



/oyi - fo Ao 



fo,x 



i^o,x - foy - fo - Jq- 

Because fo ranges from /oo to and foo > \/2/3, we 
know that fo either starts at or passes through y^2/3. 
(Strictly speaking we should be writing here /oo.Oi the 
lowest order term in t he exp ansion of /oo-) Thus, the 
effect of the pole in Eq. ( 4.24a ) must be considered. If it is 
not canceled by a zero in fio, fo,x diverges at fo = ■\/2/3. 

W e can obtain a n exp ression for Ao(/o) by dividing 
Eq. ([4.24h|) by Eq. (H.24a|), which leads to 




8 



Ao dj^o 



/oVl-/o'-^* 



2 - 3/0^ 



d/o. (4.25) 



Integrating both sides and recalling the boundary condi- 
tion = 0, we find 

_ ?2 _ 3 "4 ^ 2 1 3 „4 
2 ~ JO ^ iO J ca + ^ ioo 



+2 J* In 



-2 J* In 



/o 



1 + 



+ 3J*VWi:, (4.26) 



where Jq = fla\n~fL- To keep /o,a; from diverging, 
we insist that fiQ^fo ~ ■\/2/3) = which can be shown 
from Eq. ( 1.26| ) to imply /oo = ■\/2/3, i.e. the small-M 
limit of the stall current is the critical depairing current. 
Note that the pole in Eq. (4.24a) and the compensating 
zero in /io(/o) occur at t he bou ndary {x — oo). 
We can rearrange Eq. ( 4.24a ) as follows 



(2 - 3/2) df 



/o(o) f\/T^n^ Mf) 



(4.27) 



Then we can substitute in Eq. ( 4.26| ) for /io(/)i numeri- 
cally integrate the resulting expression and finally invert 
it in order to calculate fo{x), the u ^ profile. Fig- 
ure 1^ includes a comparison of fo{x) and the profile of a 
small-M numerical solution. 

To find the asymptotic behavior of /o and fla in the 
superconducting region, Taylor expand /to(/o) around foe 



fJ-o 



(/o) = -3V2| 



(4.28) 



Notice that /io(/o) is a second order zero, so that /o.j = 0, 
as it should at the boundary. As a conseq uenc e, the 
integral supplying the inverse profile, Eq. ( 4.27 ), has 
a pole; integrating the expression in its neighborhood 
yields \/61n(-\/2/3 — fa), leading to 



/o(i) 



Ao exp 



(4.29) 



where Aq is an integration constant undetermined be- 
cause of the translational invariance. Note /o(a;) has the 
form assumed in the preliminary analysis with A^^o = 



'u/6. Putting this result into Eq. (4.24a) leads to 



Mo 



(x) 3V2AI exp (2x/Vtj , 



(4.30) 



where A^,o = u^^^foo in agreement with the expression 
found previously. 



Let us examine Eqs. (4.24a) and ( 4.24b| ), which are 
strictly speaking superconducting solutions, in the nor- 
mal (small- /o) limit. Eq. (4.241:) leads to iJ,o{x) w —JcX, 
and inserting this into Eq. (4.24a) reveals that /o ^ in 
the following way 



fo{x) exp I 



(4.31) 



This same dependence was seen earlier in t he a nalysis of 
the bump shapes in the small-J limit, Eq. (3.6). 

What is surprising here is that what are ostensibly 
the "outer" equations for the superconducting region also 
satisfy the boundary conditions in the normal region and 
interpolate in between. This is consistent with the nu- 
merical observation that there does not seem to be a 
boundary layer at small u, that the ufxx term is appar- 
ently not a singular perturbation. With this in mind, we 
pursue the perturbative expansion to higher orders. 

The Oiu) equations. The eigenvalue Jg was deter- 
mined by examining the behavior deep in the supercon- 
ducting region and did not require imposing the bound- 
ary conditions on the normal side. Furthermore, the spa- 
tial dependence of the solution in this region is of the 
form assumed in Eqs. (4.6). We exploit these features to 
obtain higher order terms. The 0{u) equations are 

/o,x£ - IqaMi - qlh - Sfoh = 0, (4.32a) 
Ao/i + /oAi = 2qo/i,i + 2/o,sqi + foqi.x + go^A (4-32b) 
Jl = 2foqofl + f^qi - (4.32c) 

The asymptotic form of /o(x) is 

Mx) ^ /r + f^'^ G-/^ + c^-Zv^ + . . . (4.33) 



and similarly for qo{x) and fio{x). Eqs. (4.32) can be 
satisfied if the asymptotic form of fi (x) is 



(4.34) 



and similarly for qi{x) and fi,i{x). At 0{u^), f2{x) would 
have second-order polynomials multiplying the exponen- 
tials, and so on. Substituting these expressions into the 
differential equations allows us to determine the unknown 
constants (except for those associated with the transla- 
tional invariance). For /oo it yields the series 



(4.35) 



/oo \/ 2 



24 V6 

which corresponds to 
2 

J* = 



768V6 



3n/3 576V3 5184%/3 



(4.36) 



Note that the first correction to the u ^ limit of J* is 
of 0{v?), since the lowest term Jc is at the maximum of 

J*{Ioo)^flV^^. 
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The series found through the asymptotic perturbative 
expansion above can b e obt ained by a nother method. 
Looking back at Eqs. (4.29) and ( 4.3Cl| ), we note that 
the ratio of decay lengths A//A^ — 2. If we insert the 



expr essio ns we have for these length scales, Eqs. (4.7c) 
and (L9), we find as u ^ 



hi 

Xn 



ufl 



1/2 



(4.37) 



[6/^ - 4_ 

Solving for foe , and then calculating J* , we find 

J* = J,(l - u/8)i/2(l - u/24)-3/2^ (4.38) 



with Jc = \/ 4/27, w hich when expanded for small-u 
agrees with the series ( 4.36| ) found above. We plot the 
small-w numerical data and this expression together in 
Fig. ^. The fit is surprisingly good at small u, suggesting 
to us that the corrections to Eq. ( 4.38 ) are exponentially 
small as u ^ 0. 



0.010 



0.008 



0.006 



0.004 - 



0.002 - 



0.000 




FIG. 7. ( Jc — J) for the numerical d ata ( solid line) and for 
the result of the small-u analysis, Eq. (4.3S) (dashed line). 



V. MOVING INTERFACES 

At currents other than J* , the NS interfaces move with 
a constant velocity. For such solution s th e operator dt 
can be replaced by ~cdx, so that Eqs. (4.3) become 



-cufx^fxx-fq^ + f-f, 
u{-cq + fi)f = 2fxq + fq,,, 

J ^ Pq~ Ma:- 



(5.1a) 
(5.1b) 
(5.1c) 



While the boundary conditions on / and q remain the 
same, that on the scalar potential becomes — cq^a- 
Actually, it is more convenient to use instead jl = ^ — cq, 
which is the gauge-invariant potential in the constant- 
velocity case. 

The superconducting phase invades the normal phase 
if J < J* and vice versa if J > J* . For currents near J* , 



the interface speed is proportional to (J — J*). In this 
linear response regime, one can define a kinetic coefficient 
(which Likharevliy refers to as a "viscosity" ) 



dc 
dJ 



(5.2) 



J=J' 



Figure ^ shows the numerically determined kinetic coef- 
ficient as a function of u. For large-u, we find ry ^ v?^'^ 
for which we provide an argument below. 




FIG. 8. A log-log plot of the numerically determined ki- 
netic coefficient as a function of u (the solid line) along with 
an asymptotic fit of 0.797u^''* (the dotted line). 

Farther from the stall current, the velocities deviate 
from this linear behavior, as seen in Fig. |. The greatest 
departure-pccurs in the limits J — > and J — > Jc. In fact, 
LikharevtJ conjectured that the interface speed diverges 
in both of these limits; we find that it is bounded. 




0.1 0.2 0.3 

Applied current, J 

FIG. 9. The velocity of the front versus the current J for 
u = 5.79. 

The J — > limit. The moving interface equations. 
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Eqs. ( |5.l| ), simplify in the J ^ 
implies that both q ^ and /i 



limit, since that limit 
-> 0, leaving only 



fx. + ucf., + / - f = 0. 



(5.3) 



If we replace uc in the above equation by a speed w, then 
we have the steady-state version of Fisher's equationnZl, 
which is. known to have propagating front solutions with 



In our case this implies that as J 



0, 



2/u, which is in good agreement with the numerical data 
shown in Fig. ||. 

We can combine the above result with an earlier one 
to suggest that r] ~ it'^/^ as m ^ oo. In the large-w hmit, 
we have information on the following two points: (1) the 
stalled interface (J = J* ~ ii~^/^,c — 0); and (2) the 
interface in the absence of current [J — Q,c— 2/u). In 
going from (1) to (2), the changes in current and velocity 
are A J ~ u~^^^ and Ac ~ u~^. As u — > oo, both of 
these changes are small so that i] might be approximated 
by 



AJ 
A^^ 



,3/4 



(5.4) 



yielding the behavior seen in the numerical data (see 
Fig. I and Table H). 

The J ^ Jc limit. The numerical work indicates that 
the velocity is finite as J —> Jc; the limiting velocity 
is shown in Fig. ^ as a function u. We can find an 
analytic bound on this velocity as follows. First, take 
Eqs. (5.1), use the gauge-invariant p oten tial jl, and find 
the constant- velocity analog of Eq. (45). Then substi- 
tute the asymptotic forms, Eqs. ([1.6|), into the resulting 
equations, leading to 



(^/^-A;^)Ale^/^"=cglA-VA^ 
-/iiA-^e^/-^^ = 0. 



(5.5a) 
(5.5b) 



(5.5c) 



Arguments similar to those following Eqs. (4.7) lead one 



to the conclusion that in this case A/ = Ag = A^. The 
above equations can then be shown to yield the following 
relation 



2 2 
u c 



(2u\- 



2ufl\-u'fl\)c 
+ [2 {uflX^ - 1) (3/^, - 2) - ufl + A-2] = 0, (5.6) 



where we have used = 1 — /^. We find the bound by: 
(1) solving Eq. (5.6) for c; (2) substituting in /oo = \/2/3 
(which corresponds to J = Jc); and (3) extremizing that 
result with respect to the decay length A. The small-w 
limit of the resulting bound is — -s/^u/O, and the large- 
u limit is — 1/2\/3. The square-root dependence of the 
velocity in the small-w limit agrees with the data. Now 
we can consider going from the stall current (J* , c ~ 
0) to the critical depairing current ( Jc, c ^ m^/^) which 



results in changes AJ ~ v? and Ac 



v}^"^, suggesting 

that the small-w kinetic coefficient 77 ~ u^/^, which is in 
rough agreement with the numerical data. We have also 
observed that as a function of J the speed appears to 
approach its bound via a square root dependence 



c( J) = A + S( Jc - J) 



1/2 



(5.7) 



= for all 



0.20 




FIG. 10. The velocity as J —» Jc as a function of u. For 
large u, the velocity asymptotically approaches 0.9214"" *^. 



VI. SUMMARY AND REMARKS 

In this paper we have studied in detail the nucleation 
and growth of the superconducting phase in the presence 
of a current. The finite amplitude critical nuclei grow as 
the current is increased, with the amplitude eventually 
saturating as the stall current J* is approached, leading 
to the formation of interfaces separating the normal and 
superconducting phases. The stall current can be calcu- 
lated in the limit of large u using matched asymptotic 
expansions, demonstrating once again the utility of this 
technique for problems in inhomogeneous superconduc- 
tivity. We have also derived an analytic expression for 
the stall current for small m, which we believe to be cor- 
rect up to exponentially small corrections. Deviations 
from the stall current cause the interfaces to move, and 
we have calculated the mobility of these moving inter- 
faces for a wide range of u. Finally we have shown that 
the interface velocity c = 2/MasJ— >0 and that c is 
bounded as J —> Jc, in contrast to some conjectures in 
the literature. 

As in the magnetic-field analogy, the issue of stability 
and dynamics of the current-induced NS interfaces will be 
more complicated and interesting in the two-dimensional 
case. Some preliminary wack in this direction has been 
reported by Aranson et a/£3, who find that the current 
has a stabilizing effect on the NS interface. This can be 
interpreted as a positive surface tension for the interface. 
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due entirely to nonequilihrium effects. Tfiey provide a 
heuristic derivation of an interesting free-boundary prob- 
lem for the interfacial dynamics (a variant of the Lapla- 
cian growth problem); however, this free-boundary prob- 
lem is sufficiently complicated that they were unable to 
solve it to compare with their numerical results. Clearly, 
further work in this direction would be helpful in under- 
standing the nucleation and growth of the superconduct- 
ing phase in two-dimensional superconducting films. 
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APPENDIX A: AMPLITUDE OF THE CRITICAL 
NUCLEI IN THE J -> LIMIT 

In this appendix we provide a self-consistent calcula- 
tion of the amplitude of the critical nuclei in the J ^ 
limit. Choosing the gauge appro pria te for bumps cen- 
tered at X = and combining Eqs. ( |2.5| ) into one equation 
yields 

[-udt + iuJx + dl + l] il^{x,t) = \'\}j{x,t)\'^ i'{x,t) 



+iu 



dy Im {ip*{y,t)dyip{y,t)) 



^{x,t). (Al) 



The propagator for the linear operator appearing on the 
left hand side of Eq. (Al) satisfies the condition 



[-udt + iuJx + 9^ + 1] G{x, x';t- t')Q{t - t') 
= -u 5{x - x')S{t - t') 



(A2) 



and is given by 



G{x, x'; t) = 



1/2 



exp 



12m 



X exp 



ijT{x + x') u{x — x') 



M21 



4t 



(A3) 



Ivlev et aZ.L^lO used this linear propagator to evolve per- 
turbations having widths of 0(1) and carrying no cur- 
rent. Without the nonlinear terms such perturbations 
initially grow but ultimately reach a maximum size and 
then decay away. Ivlev et al. suggested that the ampli- 
tudes of the critical nuclei are exponentially small in the 
J —f limit by asking what sized initial perturbations 
are of 0(1) at their maxima. Their arguments motivated 
us to use the propagator in a more careful estimate of 
the amplitude that includes the nonlinear terms as an 
essential ingredient. 

We can convert Eq. ( |Al[ ) into an integral equation by 
multiplying both sides of Eq. ( |A1| ) (with x —i- x' and 
t t') by G{x, x' \ t — t') and integrating over all x' and 



integrating t' from to t. After some manipulation these 
steps lead to 



V'(x,i)= / dt' 
Jo 



dx' Gix,x',t~t') 



X <^ ij{x', t')5{t - t') - -\^{x', t')\^ ^{x\ t') 



dylm{r{y,t')dyij{y,t')) 



^ix',t')\, (A4) 



where t > 0. 

In order to estimate the amplitude of the threshold 
solutions, we will substitute into Eq. (A4) the following 
form 



ip{x,t) = ^/-oexp 



iJx'^ 



(A5) 



Note that this form is stationary and has a fixed Gaussian 
shap e (w hich is inspired by our WKB approximation, see 
Eq. ( |3.6D ) but it has an arbitrary amplitude which we will 
determine self-consistently. 

Let us take the t ^ oo limit and focus on x ~ 
since our interest is in the amplitude. After substitut- 
ing Eq. ( |A5| ) into Eq. (A4), we can do both integrals for 
the first term on the right hand side exactly, and it can 
be seen to decay to zero in the t —^ oo limit. Next, we 
perform the x' integration of the second term on the right 
hand side (//), which yields 



uJ JO 



dr 



exp 



24t2 - 4t3 - 3r'^ 
12uJ(1 + 3t) 



(A6) 



where t — Jt. We now apply the method of steepest 
descent to obtain 



11 exp 



32 



81 uJ 



(A7) 



In the third term on the right hand side (///) of 
Eq. (A4), wc make the substitution y = vx' and then 
perform the x' integration giving 



/// = 



j4_ 



X exp 



Jo ^[l + Til + 2v^)]^ 
24^2^2 _ _ (1 + 2v^)t^ 
12uJ[H- (1 + 2w2)r] 



(A8) 



The maximum of the term in the exponential of /// oc- 
curs at u = 1 (which is an endpoint). Linearizing about 
that maximum provides 



/// 



dTT{2 + t) r 24t2 - 4t3 - 3r4 
'^^^ \ 12uJ(l + 3r) 



uJ^ Jo V(1 + 3t)3 

X dw exp < -TT 

Jo I ^^J(l + 3r)2 



(A9) 
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where w = 1 — v. After the w integration, we apply the 
method of steepest descent to the r integration to obtain 



///: 



2J 



■ exp 



32 



81 uJ 



(AlO) 



Putting all of these results back into Eq. (A4) gives 



■00 



TTUtpQ 



■ exp 



32 



81 uJ 



(All) 



which provides the expression given in the text, Eq. ( 3 



This calculation clearly runs into trouble when u < 8/9; 
however, the numerical coefficients in front of these in- 
tegrals are sub-leading terms, and they can be varied by 
adding sub-leading terms to the initial Gaussian guess. 
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